library(tidyverse)library(tidylog)library(patchwork)library(viridis)library(RColorBrewer)library(ggforce)library(readxl)library(ggalluvial)#library(ggspatial)library(raster)library(sp)# devtools::install_github("seananderson/ggsidekick") # not on CRAN library(ggsidekick); theme_set(theme_sleek(base_size =10))# Set make working directory relative to your path -- no hardcoded pathshome <- here::here()# Load functions for plotting and manipulating datasource(paste0(home, "/R/lat-lon.R"))source(paste0(home, "/R/map-plot.R"))# Change the limit for scientific notationoptions(scipen =10000)
Load data
Making figures in the order they appear in the previous atlas
Sometimes I have standardized omitted plots and I have focused more on producing a script that can be updated each year automatically
Trends over time (single level): geom_line(), and also when the trend is not visible with bar or area plot
Trends over time (many levels): geom_area() or geom_bar, depending on how smooth it is over time. If smooth, area, if not bar. Bar also for discards because we have gaps in the time series
many levels (time periods or single time): geom_col() with facet_zoom() to visualize small groups
Seasonal trends: coord_polar()
Colors:
viridis plasma for map plots
gears: ramp up RdYlBu (cluster by fishery)
species: ramp up set3 (cluster by family)
fisheries: mako (cluster by similarity)
years: spectral
quarto, renv::renv(), GitHub to ensure reproducibility
TODO
Ansträngning? Days-at-sea? Over time or within a year?
Rerun with high-quality maps
Gör även fig. 2.1.10 och 2.1.11?
Långsiktigt: - Gör funktioner för data prep! Typ, get_vessel_size_distribution() - Värden på arter
Code
# Sofia has now updated it so that it's a single file with the "correct" focat codeload(paste0(home, "/data/test_lovgren_foCatFA.RData"))lb <- lovgren_foCatFA %>%rename(foCatFA = foCatFA_update)# Add fisheries common namesfishery_names <-read_delim(paste0(home, "/data/names_fisheries.csv"), delim =";") %>% dplyr::select(-HAVSKOD, -Redskap, -foCatFA_tmp, -foCatEu4)# See how they match#setdiff(fishery_names$foCatFA, lb$foCatFA)lb <- lb %>%left_join(fishery_names, by ="foCatFA")# TODO: All active non-pelagic gear in the Baltic are not cod-trawl, esp after 2019. Rename?# Add coordinateslb <- lb %>%mutate(lat =format_position(LATITUD),lon =format_position(LONGITUD))# Now that we have coords, we can trim the map plot function# High resolution maps for single plots and/or local plotsplot_map_l <- plot_map_l +xlim(min(filter(lb, lat >50)$lon), max(filter(lb, lat >50)$lon)) +ylim(min(filter(lb, lat >50)$lat), max(filter(lb, lat >50)$lat))# Medium resolution map is for facets and/or for large mapsplot_map_m <- plot_map_m +xlim(min(lb$lon), max(lb$lon)) +ylim(min(lb$lat), max(lb$lat))# Add year & monthlb <- lb %>%mutate(year =year(as.Date(as.character(AVRESDAT), format ="%Y%m%d")),month =month(as.Date(as.character(AVRESDAT), format ="%Y%m%d")),month_name =as.factor(month),month_name =fct_recode(month_name,"Januari"="1","Februari"="2","Mars"="3","April"="4","Maj"="5","Juni"="6","Juli"="7","Augusti"="8","September"="9","Oktober"="10","November"="11","December"="12"))# Add harbour names. For some reason they are in two files... harbors <-read_excel(paste0(home, "/data/Code-location-v1.7.xls"), skip =4) %>% dplyr::select(`ISO-3 Country Code`, `LOCODE Code`, Description, Latitude, Longitude) %>% janitor::clean_names() %>%rename(LANDHAMN_EUKOD = locode_code)# Join harbour file with coordinates to the landhamn data below (LANDHAMN_EUKOD = LOCODE, SE 999 is NA)load(paste0(home, "/data/landhamn_2003-2022.RData")) # The way to join harbours into the lb data is not via harbour name, but unique activity. That's we we join harbors to landhamnlandhamn <- landhamn %>%left_join(harbors, by ="LANDHAMN_EUKOD")# Resnummer, loggbladsnummer, aktivitetsnummer, year is an ID for a single "effort". But it's duplicated, but not in any of the columns I can see. Perhaps landed species? Make this ID in both and join. landhamn <- landhamn %>%mutate(id =paste(RESNR, LOGGBLNR, AKTIVNR, year, sep ="_")) %>%distinct(id, .keep_all =TRUE) %>%rename(landhamn_lat = latitude,landhamn_lon = longitude) %>% dplyr::select(id, LANDHAMN, landhamn_lat, landhamn_lon)lb <- lb %>%mutate(id =paste(RESNR, LOGGBLNR, AKTIVNR, year, sep ="_")) lb <- lb %>%left_join(landhamn, by ="id")# Efter 2016 finns det olika typer av landnings (för landhamn så blir det ingen match om inte fångsten landats! se t.ex. redskap)# Add species common and latin namesspcode <-read_delim(paste0(home, "/data/names_species.csv"), delim =";") %>% dplyr::select(FNAMNSVE, MAFKOD)# Clean some errors and categories we don't uselb <- lb %>%left_join(spcode, by ="MAFKOD") %>%filter(KVANT >0) %>%filter(!FNAMNSVE %in%c("Redskap Med Sälskada (St Redskap)","Obest. benfiskar")) %>%drop_na(FNAMNSVE)# Add gear names# TODO: use Sofia's data? But 2021 does not contain all gearsgear_names <-read_delim(paste0(home, "/data/scrap_names_gear.csv"), delim =";") %>% dplyr::select(-Note)# gear_names <- read_delim(paste0(home, "/data/names_gear.csv"), delim = ";") %>%# dplyr::select(-Note) %>% # mutate(Redskap2 = NA) %>% # mutate(Redskap2 = ifelse(grepl("Parbottentrål", REDSKAP), "Parbottentrål", Redskap2),# Redskap2 = ifelse(grepl("Bottentrål", REDSKAP), "Bottentrål", Redskap2),# Redskap2 = ifelse(grepl("Btrål", REDSKAP), "Bottentrål", Redskap2),# Redskap2 = ifelse(grepl("Flyttrål", REDSKAP), "Flyttrål", Redskap2),# Redskap2 = ifelse(grepl("Parflyttrål", REDSKAP), "Parflyttrål", Redskap2)) %>% # as.data.frame()# gear_names# # unique(gear_names$Redskap2)lb <- lb %>%left_join(gear_names, by ="REDSKKOD")# Load discardsdiscard <-read_excel(paste0(home, "/data/Discard 2008-2022.xlsx")) %>% dplyr::select(Year, foCatNat, Art, `Utkast (kg)`) %>%rename(year = Year,discards =`Utkast (kg)`)# Join the total loggbook landings by the same aggregation (year, fishery)lb_fishery <- lb %>%rename(Art = FNAMNSVE) %>%summarise(KVANT =sum(KVANT), .by =c(foCatNat, year, Art))discard <- discard %>%left_join(lb_fishery, by =c("foCatNat", "year", "Art")) %>%filter(discards >0)# Lastly, join in fishery based on focatnatdiscard <- discard %>%left_join(lb %>%distinct(foCatNat, .keep_all =TRUE) %>% dplyr::select("fishery", "foCatNat", "HAVSKOD"),by ="foCatNat") %>%drop_na(fishery)# Annotation df for discardmd <-expand_grid(year =c(2020, 2021, 2015), y =0,fishery =unique(discard$fishery),label ="No discard data collection") %>%mutate(keep =ifelse(year ==2015&!fishery =="Bottentrålfiske efter torsk i Östersjön", "N", "Y")) %>%filter(keep =="Y") %>% dplyr::select(-keep)# Add ICES subdivisionshape <-shapefile(paste0(home, "/data/ICES_StatRec_mapto_ICES_Areas/StatRec_map_Areas_Full_20170124.shp"))pts <-SpatialPoints(cbind(lb$lon, lb$lat), proj4string =CRS(proj4string(shape)))lb$subdiv <-over(pts, shape)$Area_27
# TODO: Check Redskap2, did I do it correctly? I think I could also have gotten redskap by joining Fiskeridata_FLEET by "foCatEu4", but then I'd rather make a new file that is distinct with respect to "foCatEu4" and "redskap", because it's unclear to me why it isn't (has year/month but they don't change by that...). Moreover, it seems to miss things. Perhaps it's best to match and get Redskap and then assign that to a category
Entire fleet
Map-plots of fishing activity (landings)
Code
# For creating hexagon polygons: https://urbandatapalette.com/post/2021-08-tessellation-sf/# For summing landings by hexagon polygons https://stackoverflow.com/questions/69438693/tallying-variable-for-points-across-polygons-with-sfst-coverslb_sf <- lb %>%st_as_sf(coords =c("lon", "lat"), crs =4326, remove =FALSE)area_honeycomb_grid =st_make_grid(lb_sf, 0.5, what ="polygons", square =FALSE)#area_square_grid = st_make_grid(lb_sf, c(0.15, 0.1), what = "polygons", square = FALSE)# To sf and add grid IDhoneycomb_grid_sf =st_sf(area_honeycomb_grid) %>%mutate(polygon_id =as.character(1:length(lengths(area_honeycomb_grid)))) # add grid IDpolygon_counts <- lb_sf %>%mutate(polygon_id =as.character(st_intersects(geometry, honeycomb_grid_sf))) %>%group_by(year, polygon_id) %>%summarise(sum_kvant =sum(KVANT)) %>%st_drop_geometry()polygon_sum <- honeycomb_grid_sf %>%left_join(polygon_counts, by ="polygon_id") %>%drop_na(sum_kvant)# All years but Morocco & Arctic excluded for visibilitylandings_map <- plot_map_m +geom_sf(data = polygon_sum, aes(fill = sum_kvant/1000), color =NA) +scale_fill_viridis(trans ="sqrt", option ="plasma",na.value ="yellow",limits =c(0, quantile(polygon_sum$sum_kvant/1000, 0.99)),name ="Landningar (ton)") +geom_sf(color ="gray70") +facet_wrap(~year, ncol =5) +xlim(min(filter(lb, lat >50)$lon), max(filter(lb, lat >50)$lon)) +ylim(min(filter(lb, lat >50)$lat), 70) +guides(fill =guide_colorbar(title.position ="top", title.hjust =0.5)) +theme(axis.text.x =element_text(angle =90),legend.position ="bottom", legend.title =element_text(size =6),legend.text =element_text(size =5),legend.key.width =unit(1, "cm"),legend.key.height =unit(0.3, "cm"))ggsave(plot = landings_map, paste0(home, "/figures/entire_fleet/landings_map.png"), width =17, height =17, units ="cm")# Now aggregate all years pooled and make a big map plotpolygon_counts_agg <- polygon_counts %>%mutate(tp =ifelse(year >2015, "2016-2022", "2003-2015")) %>%group_by(tp, polygon_id) %>%summarise(sum_kvant =sum(sum_kvant) /1000000)polygon_sum_agg <- honeycomb_grid_sf %>%left_join(polygon_counts_agg, by ="polygon_id") %>%drop_na(sum_kvant)landings_map_select <- plot_map_m +geom_sf(data = polygon_sum_agg, aes(fill = sum_kvant), color =NA) +scale_fill_viridis(trans ="sqrt", option ="plasma",na.value ="yellow", name ="Landningar (tusen ton)",limits =c(0, quantile(polygon_sum_agg$sum_kvant, 0.99))) +geom_sf(color ="gray70") +facet_wrap(~tp) +theme(axis.text.x =element_text(angle =0),legend.position =c(0.12, 0.9))ggsave(plot = landings_map_select, paste0(home, "/figures/entire_fleet/landings_map_select.png"), width =17, height =17, units ="cm")
Fig. x Havsområde, art, landningsland. Topp 10 arter, och för dessa, topp 7 landningsländer
Cumulative number of harbours
Code
n_harbours <- lb %>%summarise(n_harb =length(unique(LANDHAMN)),n_vessel =length(unique(BATNAMN)), .by = year) %>%mutate(n_harb_scaled = n_harb / n_vessel)p1 <-ggplot(n_harbours, aes(year, n_harb)) +geom_line(color ="gray40") +labs(x ="År", y ="Antal landningshamnar")p2 <-ggplot(n_harbours, aes(year, n_harb_scaled)) +geom_line(color ="gray40") +labs(x ="År", y ="Antal landningshamnar relativt antal fartyg")harbors <- p1 / p2 +plot_layout(axes ="collect")ggsave(plot = harbors, "figures/entire_fleet/harbors.png", width =17, height =15, units ="cm")cumsum_harb <- lb %>%summarise(kvant_sum =sum(KVANT),.by =c(year, LANDHAMN)) %>%arrange(year, desc(kvant_sum)) %>%mutate(cum_sum =cumsum(kvant_sum), year_tot =sum(kvant_sum),cum_sum_prop = cum_sum / year_tot,n =row_number(),.by = year)ogives_harb <- cumsum_harb %>%filter(year %in%c(2003, 2022) & cum_sum_prop <0.8) %>%group_by(year) %>%filter(cum_sum_prop ==max(cum_sum_prop))# Here I need a color that has visible colors in both extreme endsnyears <-length(unique(cumsum_harb$year))colors <-colorRampPalette(brewer.pal(name ="Spectral", n =10))(nyears)p1 <-ggplot(cumsum_harb, aes(n, cum_sum_prop, color =factor(year))) +geom_line(alpha =0.9) +scale_color_manual(values = colors, name ="") +labs(y ="Kumulativ andel", x ="") +theme_sleek(base_size =9) +theme(legend.position =c(0.65, 0.1),legend.direction ="horizontal",legend.spacing.y =unit(0, 'cm'),legend.spacing.x =unit(0, 'cm')) +guides(color =guide_legend(nrow =2))p2 <-ggplot(cumsum_harb, aes(n, cum_sum_prop, color =factor(year))) +geom_hline(yintercept =0.8, alpha =0.5, linetype =2) +geom_line(alpha =0.9) +xlim(0, 65) +geom_segment(data = ogives_harb, aes(x = n, xend = n, y =0, yend = cum_sum_prop, color =factor(year)),alpha =0.9, linetype =2) +scale_color_manual(values = colors, name ="") +labs(y ="Kumulativ andel", x ="Antal landningshamnar") +theme_sleek(base_size =9) +theme(legend.position ="bottom") +guides(color ="none") +coord_cartesian(expand =0)cumulative_harbour <- p1 / p2ggsave(plot = cumulative_harbour, "figures/entire_fleet/cumulative_harbour.png", width =17, height =15, units ="cm")
Fig. x Antal landningshamnar (totalt och relativt till antal fartyg)
Fig. x Kumulativ andel av landningsvikt för det svenska yrkesfisket per år och landningshamn. De vertikala linjerna i den nedre panelen visar det antal hamnar som tar emot 80% av landningsvikten för åren 2003 (röd) och 2022 (blå).
Fig. 1.4.2–1.4.3 Längdfördelning av fartyg den svenska fiskeflottan per år (överst), förändring över tid i andel och procent i fartyg per storleksklass (nedre vänster) samt förändring över tid totalt över alla storlekar, över och under 10m (nedre höder)
Fig. 1.4.12 Antal fartyg per fiske i det svenska yrkesfisket på kusten och i havet för 2003-2015 (ljus) och 2016-2022 (mörk). Pilar indikerar skillnad och riktning
Cumulative share of landings
Code
cumsum <- lb %>%summarise(kvant_sum =sum(KVANT),.by =c(year, BATNAMN)) %>%arrange(year, desc(kvant_sum)) %>%mutate(cum_sum =cumsum(kvant_sum), year_tot =sum(kvant_sum),cum_sum_prop = cum_sum / year_tot,n =row_number(),.by = year)ogives <- cumsum %>%filter(year %in%c(2003, 2022) & cum_sum_prop <0.8) %>%group_by(year) %>%filter(cum_sum_prop ==max(cum_sum_prop))# Here I need a color that has visible colors in both extreme endsnyears <-length(unique(cumsum$year))colors <-colorRampPalette(brewer.pal(name ="Spectral", n =10))(nyears)p1 <-ggplot(cumsum, aes(n, cum_sum_prop, color =factor(year))) +geom_line(alpha =0.9) +scale_color_manual(values = colors, name ="") +labs(y ="Kumulativ andel", x ="") +theme_sleek(base_size =9) +theme(legend.position =c(0.65, 0.1),legend.direction ="horizontal",legend.spacing.y =unit(0, 'cm'),legend.spacing.x =unit(0, 'cm')) +guides(color =guide_legend(nrow =2))p2 <-ggplot(cumsum, aes(n, cum_sum_prop, color =factor(year))) +geom_hline(yintercept =0.8, alpha =0.5, linetype =2) +geom_line(alpha =0.9) +xlim(0, 65) +geom_segment(data = ogives, aes(x = n, xend = n, y =0, yend = cum_sum_prop, color =factor(year)),alpha =0.9, linetype =2) +scale_color_manual(values = colors, name ="") +labs(y ="Kumulativ andel", x ="Antal fartyg") +theme_sleek(base_size =9) +theme(legend.position ="bottom") +guides(color ="none") +coord_cartesian(expand =0)cumulative_landings <- p1 / p2ggsave(plot = cumulative_landings, "figures/entire_fleet/cumulative_landings.png", width =17, height =15, units ="cm")
Fig. 1.4.4 Kumulativ andel av landningsvikt för det svenska yrkesfisket per år. De vertikala linjer i den nedre panelen visar det antal fartyg som fiskar in 80% av landningsvikten för åren 2003 (röd) och 2022 (blå).
Landings by gear
Code
lb_gear <- lb %>%drop_na(Redskap2) %>%summarise(sum_kvant =sum(KVANT), .by =c(year, Redskap2))pal_long <-colorRampPalette(brewer.pal(name ="RdYlBu", n =11))(length(unique(lb_gear$Redskap2)))order <- lb_gear %>%summarise(sum_kvant_all =sum(sum_kvant), .by = Redskap2) %>%arrange(desc(sum_kvant_all)) %>%mutate(sum =sum(sum_kvant_all),percent =100* (sum_kvant_all / sum)) %>%mutate(group =ifelse(percent >1, ">1% av totala landningarna", "<1% av totala landningarna"))lb_gear <- lb_gear %>%left_join(order, by ="Redskap2")p1 <-ggplot(lb_gear %>%filter(group ==">1% av totala landningarna"),aes(year, sum_kvant/1000000, fill = Redskap2, order = sum_kvant)) +#geom_bar(stat = "identity") +geom_area() +scale_fill_manual(values = pal_long, name ="") +theme_sleek(base_size =9) +theme(legend.position =c(0.8, 0.93),legend.direction ="horizontal",legend.spacing.y =unit(0, 'cm'),legend.spacing.x =unit(0.2, 'cm'),legend.key.size =unit(0.4, 'cm')) +labs(y ="Landningar (tusen ton)", x ="") +guides(fill =guide_legend(nrow =2)) +coord_cartesian(expand =0) +facet_wrap(~group)p2 <-ggplot(lb_gear %>%filter(group =="<1% av totala landningarna"),aes(year, sum_kvant/1000000, fill = Redskap2, order = sum_kvant)) +#geom_bar(stat = "identity") +geom_area() +scale_fill_manual(values =rev(pal_long), name ="") +theme_sleek(base_size =9) +theme(legend.position =c(0.68, 0.93),legend.direction ="horizontal",legend.spacing.y =unit(0, 'cm'),legend.spacing.x =unit(0.2, 'cm'),legend.key.size =unit(0.4, 'cm')) +labs(y ="Landningar (tusen ton)", x ="År") +guides(fill =guide_legend(nrow =2)) +coord_cartesian(expand =0) +facet_wrap(~group)landings_gear <- p1 / p2ggsave(plot = landings_gear, "figures/entire_fleet/landings_gear.png", width =17, height =20, units ="cm")
Fig. 1.4.5 Redskap som används i det svenska yrkesfisket på kusten och i havet över tid. Andelen är beräknat på landningsvikt.
Fig. 1.4.7 De vanligaste 12 landade arter i det svenska yrkesfisket på kusten och i havet över tid.
Fig. 1.4.8 De vanligaste 30 landade arterna (2016-2022) i det svenska yrkesfisket på kusten och i havet.
# TODO: By value also? all species or just the ones you can actually see?
Landings map for top 12 species
Code
# Find which are the top 5 specieslb_spec_12 <- lb %>%drop_na(FNAMNSVE) %>%summarise(sum_kvant =sum(KVANT), .by = FNAMNSVE) %>%arrange(desc(sum_kvant))polygon_counts_sp <- lb_sf %>%filter(FNAMNSVE %in%c(lb_spec_12$FNAMNSVE[1:12])) %>%filter(year >2015) %>%mutate(polygon_id =as.character(st_intersects(geometry, honeycomb_grid_sf))) %>%group_by(polygon_id, FNAMNSVE) %>%# FIXME: want to do it for another year?summarise(sum_kvant =sum(KVANT)) %>%st_drop_geometry()polygon_sum_sp <- honeycomb_grid_sf %>%left_join(polygon_counts_sp, by ="polygon_id") %>%drop_na(sum_kvant)# Add proportion by species to make the spatial pattern more visiblepolygon_sum_sp <- polygon_sum_sp %>%mutate(sum_kvant_prop = sum_kvant /max(sum_kvant), .by = FNAMNSVE)# All years but Morocco & Arctic excluded for visibilitylandings_sp_map <- plot_map_m +geom_sf(data = polygon_sum_sp, aes(fill = sum_kvant_prop), color =NA) +scale_fill_viridis(trans ="sqrt", option ="plasma",na.value ="yellow",limits =c(0, quantile(polygon_sum_sp$sum_kvant_prop, 0.99)) ) +geom_sf(color ="gray70") +facet_wrap(~FNAMNSVE) +# TODO other way to trim plot?xlim(-8.2, 25.5) +ylim(51.5, 68) +guides(fill ="none") +theme(axis.text.x =element_text(angle =90))ggsave(plot = landings_sp_map, paste0(home, "/figures/entire_fleet/landings_sp_map.png"), width =17, height =17, units ="cm")
Fig. 1.4.9 Utbredning av fiskeplatser (rapporterade landningar) för de 12 ekonomiskt viktigaste arterna i svenska yrkesfisket på kusten och i havet 2019-2022. Skalan är relativ och inte jämförbar mellan arter i absoluta tal
# TODO: Same as above but for top 12 commercially?